Take-home Exercise 01

Published

February 25, 2024

1.0 Background

Human mobility, the movement of human beings in space and time, reflects the spatial-temporal characteristics of human behavior. With the advancement Information and Communication Technologies (ICT) especially smart phone, a large volume of data related to human mobility have been collected. By using appropriate GIS analysis methods, these data are potentially useful in supporting smart city planning and management.

In Singapore, one of the important source of data related to human mobility is from Land Transport Authority (LTA) DataMall. Two data sets related to human mobility are provided by the portal, they are: Passenger Volume by Origin Destination Train Stations and Passenger Volume by Origin Destination Bus Stops. One of the limitation of these data sets is that their location are biased to either bus stops or MRT/LRT stations. In 2020, another very interesting human mobility data set called Grab Posisi was released by GRAB, one of the largest shared taxi operator in South-east Asia. There are two data sets been released and one of them is for Singapore.

2.0 Objectives

Geospatial analytics hold tremendous potential to address complex problems facing society. In this study, you are tasked to apply appropriate spatial point patterns analysis methods to discover the geographical and spatio-temporal distribution of Grab hailing services locations in Singapore.

3.0 The Task

The specific tasks of this take-home exercise are as follows:

  • Using appropriate function of sf and tidyverse, preparing the following geospatial data layer in sf tibble data.frames:

    • Grab taxi location points either by origins or destinations.

    • Road layer within Singapore excluding outer islands.

    • Singapore boundary layer excluding outer islands

  • Using the extracted data, derive traditional Kernel Density Estimation layers.

  • Using the extracted data, derive either Network Kernel Density Estimation (NKDE) or Temporal Network Kernel Density Estimation (TNKDE)

  • Using appropriate tmap functions, display the kernel density layers on openstreetmap of Singapore.

  • Describe the spatial patterns revealed by the kernel density maps.

4.0 The Data

Apstial data

For the purpose of this assignment, Grab-Posisi of Singapore will be used.

Geospatial data

  • Road data set from OpenStreetMap of Geofabrik download server. The Malaysia, Singapore, and Brunei coverage should be downloaded.

  • Master Plan 2019 Subzone Boundary (No Sea) from Data.gov.sg.

5.0. Importing Packages

The R packages used in this project are:

  • sf, a relatively new R package specially designed to import, manage and process vector-based geospatial data in R.

  • tidyverse for performing data science tasks such as importing, wrangling and visualising data.

  • spatstat, which has a wide range of useful functions for point pattern analysis. In this hands-on exercise, it will be used to perform 1st- and 2nd-order spatial point patterns analysis and derive kernel density estimation (KDE) layer.

  • raster which reads, writes, manipulates, analyses and model of gridded spatial data (i.e. raster). In this hands-on exercise, it will be used to convert image output generate by spatstat into raster format.

  • maptools which provides a set of tools for manipulating geographic data. In this hands-on exercise, we mainly use it to convert Spatial objects into ppp format of spatstat.

  • tmap which provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API.

  • spNetwork, which provides functions to perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network. It also can be used to build spatial matrices (‘listw’ objects like in ‘spdep’ package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances.

  • rgdal, which provides bindings to the ‘Geospatial’ Data Abstraction Library (GDAL) (>= 1.11.4) and access to projection/transformation operations from the PROJ library. In this exercise, rgdal will be used to import geospatial data in R and store as sp objects.

  • sp, which provides classes and methods for dealing with spatial data in R. In this exercise, it will be used to manage SpatialPointsDataFrame and SpatiaLinesDataFrame, and for performing projection transformation.

Code
pacman::p_load(sp, sf, spNetwork, tmap, classInt, viridis, tidyverse, spatstat, sfdep, raster, maptools, arrow, lubridate)

6.0 Importing Data in R environment

6.1 Aspatial Data

6.1.1 Import GrabPosisi Data

Code
df <- read_parquet("data/Apstial/GrabPosisi/part-00000-8bbff892-97d2-4011-9961-703e38972569.c000.snappy.parquet")

6.2 Geospatial Data

6.2.1 Import Master Plan 2019 subzone Boundary (No Sea) Data

Code
mpsz_sf <- st_read(dsn= "data/Geospatial/MPSZ-2019", layer="MPSZ-2019")
Reading layer `MPSZ-2019' from data source 
  `C:\j00b00\IS415-GAA\Take-home_EX\Take-home_Ex01\data\Geospatial\MPSZ-2019' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

6.2.2 Import Road data

Code
roads <- st_read(dsn = "data/Geospatial/malaysia-singapore-brunei-latest-free", layer = "gis_osm_roads_free_1")
Reading layer `gis_osm_roads_free_1' from data source 
  `C:\j00b00\IS415-GAA\Take-home_EX\Take-home_Ex01\data\Geospatial\malaysia-singapore-brunei-latest-free' 
  using driver `ESRI Shapefile'
Simple feature collection with 1759836 features and 10 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 99.66041 ymin: 0.8021131 xmax: 119.2601 ymax: 7.514393
Geodetic CRS:  WGS 84

7.0 Data Wrangling

Before we proceed with our analysis, we need to preprocess the data to the correct format and clean it.

7.1 Assign the correct CRS

7.1.1 mpsz_sf

Code
mpsz3414 <- st_transform(mpsz_sf, crs = 3414)
st_crs(mpsz3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

7.1.2 roads

Code
roads3414 <- st_transform(roads, crs = 3414)
st_crs(roads3414)
Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

7.2 Convert data type of pingtimestamp from character to date-time

7.2.1 df

Code
df$pingtimestamp <- as_datetime(df$pingtimestamp)

7.2.2 Extract Grab taxi location points by origins

Code
origin_df <- df %>%
  group_by(trj_id) %>%
  arrange(desc(pingtimestamp)) %>%
  filter(row_number()==1) %>%
  mutate(weekday = wday(pingtimestamp, label=TRUE, abbr=TRUE), 
         start_hr = factor(hour(pingtimestamp)), 
         day = factor(mday(pingtimestamp)))

7.2.3 Creating a simple feature data frame from origin_df

Code
origin_sf <- st_as_sf(origin_df,
                      coords = c("rawlng", "rawlat"),
                      crs = 4326) %>%
  st_transform(crs = 3414)

7.3 Extract Singapore boundary layer excluding outer islands

7.3.1 First, filter out the outer islands

Code
outer_island <- mpsz_sf[mpsz_sf$PLN_AREA_N == "SOUTHERN ISLANDS" | mpsz_sf$PLN_AREA_N == "NORTH-EASTERN ISLANDS" | mpsz_sf$PLN_AREA_N == "WESTERN ISLANDS",]
Code
plot(st_geometry(outer_island))

Code
outer_island<- st_transform(outer_island, crs = 3414)

7.3.3 Filter out the outer island

Code
main_island <- st_difference(st_union(mpsz3414),st_union(outer_island))
Code
plot(st_geometry(main_island))

7.4 Filter out roads from roads3414 that have max speed value above 0.

Code
main_road <- roads3414 %>% filter(maxspeed > 0)

7.4.1 Find the roads that are only in Singapore excluding outer islands

Code
sg_roads <- st_intersection(main_road, main_island)

7.5 Road layer

Code
tm_shape(main_island) + tm_polygons() + tm_shape(sg_roads) + tm_lines(col="red", siz = 0.03)

7.6 Grab origin layer

Code
tm_shape(origin_sf) + tm_dots(siz = 0.03)

7.7 Converting sf data frames to sp’s Spatial* class

spatstat requires the analytical data in ppp object form. There is no direct way to convert a Spatial* classes into ppp object. We need to convert the Spatial classes* into Spatial object first.

Code
origin <- as_Spatial(origin_sf)
mpsz <- as_Spatial(mpsz3414)
main_island <- as_Spatial(main_island)
sg_roads <- as_Spatial(sg_roads)

7.8 Converting Spatial* class into generic sp format

Code
origin_sp <- as(origin, "SpatialPoints")
main_island_sp <- as(main_island, "SpatialPolygons")
sg_roads_sp <- as(sg_roads, "SpatialLines")
origin_sp
class       : SpatialPoints 
features    : 28000 
extent      : 3638.685, 50024.92, 25350.05, 49469.41  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 

7.9 Converting generic sp format into spatstat’s ppp format

convert the spatial data into spatstat’s ppp object format.

Code
origin_ppp <- as(origin_sp, "ppp")

7.9.1 Plot origin_ppp for visualisation

Code
plot(origin_ppp)

Code
summary(origin_ppp)
Planar point pattern:  28000 points
Average intensity 2.502667e-05 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Window: rectangle = [3638.69, 50024.92] x [25350.05, 49469.41] units
                    (46390 x 24120 units)
Window area = 1118810000 square units

7.10 Handling Duplicate Point

check the duplication in a ppp object

Code
any(duplicated(origin_ppp))
[1] FALSE

7.11 Creating owin object

When analysing spatial point patterns, it is a good practice to confine the analysis with a geographical area like Singapore boundary. In spatstat, an object called owin is specially designed to represent this polygonal region.

Code
main_island_owin <- as(main_island, "owin")

7.11.1 Display main_island_owin object

Code
plot(main_island_owin)

Code
summary(main_island_owin)
Window: polygonal boundary
37 separate polygons (29 holes)
                  vertices         area relative.area
polygon 1               71  5.63061e+03      8.47e-06
polygon 2               10  1.99717e+02      3.01e-07
polygon 3            12667  6.63014e+08      9.98e-01
polygon 4 (hole)         3 -3.41897e-05     -5.14e-14
polygon 5 (hole)        23 -1.99656e+01     -3.00e-08
polygon 6 (hole)        35 -1.38385e+02     -2.08e-07
polygon 7 (hole)        19 -4.39650e+00     -6.62e-09
polygon 8 (hole)       270 -1.21455e+03     -1.83e-06
polygon 9 (hole)         3 -4.95057e-02     -7.45e-11
polygon 10 (hole)        3 -3.65499e-03     -5.50e-12
polygon 11 (hole)       38 -7.79904e+03     -1.17e-05
polygon 12 (hole)        3 -5.99535e-04     -9.02e-13
polygon 13 (hole)        3 -3.04561e-04     -4.58e-13
polygon 14 (hole)        3 -7.43616e-06     -1.12e-14
polygon 15 (hole)        6 -8.37554e-01     -1.26e-09
polygon 16 (hole)        4 -2.86396e-01     -4.31e-10
polygon 17 (hole)        3 -1.81439e-04     -2.73e-13
polygon 18 (hole)        3 -8.68789e-04     -1.31e-12
polygon 19 (hole)        3 -4.46076e-04     -6.71e-13
polygon 20 (hole)        3 -3.39794e-04     -5.11e-13
polygon 21 (hole)      317 -5.11280e+04     -7.69e-05
polygon 22 (hole)        5 -2.92235e-04     -4.40e-13
polygon 23 (hole)        3 -4.52043e-05     -6.80e-14
polygon 24 (hole)        3 -3.90173e-05     -5.87e-14
polygon 25 (hole)        5 -2.44411e-04     -3.68e-13
polygon 26 (hole)        4 -2.18616e-04     -3.29e-13
polygon 27 (hole)        4 -4.28453e-01     -6.45e-10
polygon 28 (hole)        4 -2.54488e-04     -3.83e-13
polygon 29 (hole)        3 -9.59850e-05     -1.44e-13
polygon 30 (hole)       41 -4.01660e+04     -6.04e-05
polygon 31 (hole)        3 -4.14099e-04     -6.23e-13
polygon 32 (hole)        5 -3.64686e-02     -5.49e-11
polygon 33              30  2.80002e+04      4.21e-05
polygon 34              27  1.50315e+04      2.26e-05
polygon 35             285  1.61128e+06      2.42e-03
polygon 36              91  1.49663e+04      2.25e-05
polygon 37              71  8.18750e+03      1.23e-05
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units
                     (53270 x 28810 units)
Window area = 664597000 square units
Fraction of frame area: 0.433

7.12 Combining point events object and owin object

Extract Grab origin point that are in main island in Singapore

Code
originSG_ppp = origin_ppp[main_island_owin]
Code
summary(originSG_ppp)
Planar point pattern:  27821 points
Average intensity 4.186147e-05 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Window: polygonal boundary
37 separate polygons (29 holes)
                  vertices         area relative.area
polygon 1               71  5.63061e+03      8.47e-06
polygon 2               10  1.99717e+02      3.01e-07
polygon 3            12667  6.63014e+08      9.98e-01
polygon 4 (hole)         3 -3.41897e-05     -5.14e-14
polygon 5 (hole)        23 -1.99656e+01     -3.00e-08
polygon 6 (hole)        35 -1.38385e+02     -2.08e-07
polygon 7 (hole)        19 -4.39650e+00     -6.62e-09
polygon 8 (hole)       270 -1.21455e+03     -1.83e-06
polygon 9 (hole)         3 -4.95057e-02     -7.45e-11
polygon 10 (hole)        3 -3.65499e-03     -5.50e-12
polygon 11 (hole)       38 -7.79904e+03     -1.17e-05
polygon 12 (hole)        3 -5.99535e-04     -9.02e-13
polygon 13 (hole)        3 -3.04561e-04     -4.58e-13
polygon 14 (hole)        3 -7.43616e-06     -1.12e-14
polygon 15 (hole)        6 -8.37554e-01     -1.26e-09
polygon 16 (hole)        4 -2.86396e-01     -4.31e-10
polygon 17 (hole)        3 -1.81439e-04     -2.73e-13
polygon 18 (hole)        3 -8.68789e-04     -1.31e-12
polygon 19 (hole)        3 -4.46076e-04     -6.71e-13
polygon 20 (hole)        3 -3.39794e-04     -5.11e-13
polygon 21 (hole)      317 -5.11280e+04     -7.69e-05
polygon 22 (hole)        5 -2.92235e-04     -4.40e-13
polygon 23 (hole)        3 -4.52043e-05     -6.80e-14
polygon 24 (hole)        3 -3.90173e-05     -5.87e-14
polygon 25 (hole)        5 -2.44411e-04     -3.68e-13
polygon 26 (hole)        4 -2.18616e-04     -3.29e-13
polygon 27 (hole)        4 -4.28453e-01     -6.45e-10
polygon 28 (hole)        4 -2.54488e-04     -3.83e-13
polygon 29 (hole)        3 -9.59850e-05     -1.44e-13
polygon 30 (hole)       41 -4.01660e+04     -6.04e-05
polygon 31 (hole)        3 -4.14099e-04     -6.23e-13
polygon 32 (hole)        5 -3.64686e-02     -5.49e-11
polygon 33              30  2.80002e+04      4.21e-05
polygon 34              27  1.50315e+04      2.26e-05
polygon 35             285  1.61128e+06      2.42e-03
polygon 36              91  1.49663e+04      2.25e-05
polygon 37              71  8.18750e+03      1.23e-05
enclosing rectangle: [2667.54, 55941.94] x [21448.47, 50256.33] units
                     (53270 x 28810 units)
Window area = 664597000 square units
Fraction of frame area: 0.433

7.12.1 Plot originSg_ppp

Code
plot(originSG_ppp)

8.0 First-order Spatial Point Patterns Analysis

Deriving kernel density estimation (KDE) layer for visualising and exploring the intensity of point processes

8.1 Kernel Density Estimation

8.1.1 Compute kernel density estimation using automatic origin points in main land in Singapore

Code
kde_originSG_bw = density(originSG_ppp, 
                          sigma = bw.diggle, 
                          edge = TRUE, 
                          kernel="gaussian")

8.1.2 Display the kernel density derived

Code
plot(kde_originSG_bw)

8.1.3 Rescalling KDE Values

Code
originSG_ppp.km <-rescale(originSG_ppp, 1000, "km")

Re-run density() using the resale data set and plot the output kde map.

Code
kde_originSG.bw <- density(originSG_ppp.km, 
                           sigma=bw.diggle, 
                           edge=TRUE, 
                           kernel="gaussian")

8.1.4 Display the kernel density derived

Code
plot(kde_originSG.bw)

8.2 Converting KDE output into grid object Convert it so that it is suitable for mapping purposes

Code
gridded_kde_originSG_bw <- as.SpatialGridDataFrame.im(kde_originSG.bw)
spplot(gridded_kde_originSG_bw)

8.2.1 Converting gridded output into raster Convert the gridded kernal density objects into RasterLayer object

Code
kde_originSG_bw_raster <- raster(gridded_kde_originSG_bw)

Properties of kde_originSG_bw_raster RasterLayer

Code
kde_originSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4162063, 0.2250614  (x, y)
extent     : 2.667538, 55.94194, 21.44847, 50.25633  (xmin, xmax, ymin, ymax)
crs        : NA 
source     : memory
names      : v 
values     : -3.709752e-13, 1953.626  (min, max)

8.2.2 Assigning projection system in kde_originSG_bw_raster RasterLayer

Code
projection(kde_originSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_originSG_bw_raster
class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4162063, 0.2250614  (x, y)
extent     : 2.667538, 55.94194, 21.44847, 50.25633  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : v 
values     : -3.709752e-13, 1953.626  (min, max)

8.2.3 Visualising the output in tmap

Code
tm_shape(kde_originSG_bw_raster) + tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"),frame = FALSE)

8.3 Comparing Spatial Point Patterns using KDE Compare KDE of childcare at Ponggol, Tampines, Chua Chu Kang and Jurong West planning areas.

8.3.1 Extracting planning area boundary

Code
pg = mpsz[mpsz@data$PLN_AREA_N == "PUNGGOL",]
tm = mpsz[mpsz@data$PLN_AREA_N == "TAMPINES",]
ck = mpsz[mpsz@data$PLN_AREA_N == "CHOA CHU KANG",]
jw = mpsz[mpsz@data$PLN_AREA_N == "JURONG WEST",]

8.3.2 Plotting target planning areas

Code
plot(pg, main = "Ponggol")

Code
plot(tm, main = "Tampines")

Code
plot(ck, main = "Choa Chu Kang")

Code
plot(jw, main = "Jurong West")

8.3.3 Converting the spatial point data frame into generic sp format

Code
pg_sp = as(pg, "SpatialPolygons")
tm_sp = as(tm, "SpatialPolygons")
ck_sp = as(ck, "SpatialPolygons")
jw_sp = as(jw, "SpatialPolygons")

8.3.4 Creating owin object

Code
pg_owin = as(pg_sp, "owin")
tm_owin = as(tm_sp, "owin")
ck_owin = as(ck_sp, "owin")
jw_owin = as(jw_sp, "owin")

8.3.5 Combining origin points and the study area Extract childcare that is within the specific region to do our analysis later on

Code
origin_pg_ppp = originSG_ppp[pg_owin]
origin_tm_ppp = originSG_ppp[tm_owin]
origin_ck_ppp = originSG_ppp[ck_owin]
origin_jw_ppp = originSG_ppp[jw_owin]

Next, rescale() function is used to transform the unit of measurement from metre to kilometer.

Code
origin_pg_ppp.km = rescale(origin_pg_ppp, 1000, "km")
origin_tm_ppp.km = rescale(origin_tm_ppp, 1000, "km")
origin_ck_ppp.km = rescale(origin_ck_ppp, 1000, "km")
origin_jw_ppp.km = rescale(origin_jw_ppp, 1000, "km")

Plot these four study areas and the locations of the origin points

Code
plot(origin_pg_ppp.km, main="Punggol")

Code
plot(origin_tm_ppp.km, main="Tampines")

Code
plot(origin_ck_ppp.km, main="Choa Chu Kang")

Code
plot(origin_jw_ppp.km, main="Jurong West")

8.3.6 Computing KDE

Code
plot(density(origin_pg_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Punggol")

Code
plot(density(origin_tm_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Tampines")

Code
plot(density(origin_ck_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Choa Chu Kang")

Code
plot(density(origin_jw_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian"),
     main="Jurong West")

8.3.7 Storing KDE values

Code
kde_origin_pg.bw = density(origin_pg_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian")

kde_origin_tm.bw = density(origin_tm_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian")
             
kde_origin_ck.bw = density(origin_ck_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian")
             
kde_origin_jw.bw = density(origin_jw_ppp.km, 
             sigma=bw.diggle, 
             edge=TRUE, 
             kernel="gaussian")

8.4 Converting KDE output into grid object Convert it so that it is suitable for mapping purposes

Code
gridded_kde_origin_pg_bw <- as.SpatialGridDataFrame.im(kde_origin_pg.bw)
spplot(gridded_kde_origin_pg_bw)

Code
gridded_kde_origin_tm_bw <- as.SpatialGridDataFrame.im(kde_origin_tm.bw)
spplot(gridded_kde_origin_tm_bw)

Code
gridded_kde_origin_ck_bw <- as.SpatialGridDataFrame.im(kde_origin_ck.bw)
spplot(gridded_kde_origin_ck_bw)

Code
gridded_kde_origin_jw_bw <- as.SpatialGridDataFrame.im(kde_origin_jw.bw)
spplot(gridded_kde_origin_jw_bw)

8.4.1 Converting gridded output into raster

Code
kde_origin_pg_bw_raster <- raster(gridded_kde_origin_pg_bw)
kde_origin_tm_bw_raster <- raster(gridded_kde_origin_tm_bw)
kde_origin_ck_bw_raster <- raster(gridded_kde_origin_ck_bw)
kde_origin_jw_bw_raster <- raster(gridded_kde_origin_jw_bw)

8.4.2 Assigning projection system in all raster layer

Code
projection(kde_origin_pg_bw_raster) <- CRS("+init=EPSG:3414")
projection(kde_origin_tm_bw_raster) <- CRS("+init=EPSG:3414")
projection(kde_origin_ck_bw_raster) <- CRS("+init=EPSG:3414")
projection(kde_origin_jw_bw_raster) <- CRS("+init=EPSG:3414")

8.4.3 Visualizing the output in tmap

Code
tm_shape(kde_origin_pg_bw_raster) + tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"),frame = FALSE)

Code
tm_shape(kde_origin_tm_bw_raster) + tm_raster("v") + 
  tm_layout(legend.position = c("right", "bottom"),frame = FALSE)

Code
tm_shape(kde_origin_ck_bw_raster) + tm_raster("v") + 
  tm_layout(legend.position = c("right", "bottom"),frame = FALSE)

Code
tm_shape(kde_origin_jw_bw_raster) + tm_raster("v") + 
  tm_layout(legend.position = c("right", "bottom"),frame = FALSE)

9.0 Network Constrained KDE (NetKDE) Analysis Perform NetKDE analysis by using appropriate functions provided in spNetwork package.

9.1 Preparing data 9.1.1 road layer on 4 planning area

Code
pg <- mpsz3414[mpsz3414$PLN_AREA_N == "PUNGGOL",]
tm <- mpsz3414[mpsz3414$PLN_AREA_N == "TAMPINES",]
ck <- mpsz3414[mpsz3414$PLN_AREA_N == "CHOA CHU KANG",]
jw <- mpsz3414[mpsz3414$PLN_AREA_N == "JURONG WEST",]
Code
pg_roads <- st_intersection(roads3414,pg)
tm_roads <- st_intersection(roads3414,tm)
ck_roads <- st_intersection(roads3414,ck)
jw_roads <- st_intersection(roads3414,jw)

9.1.2 origin data on 4 planning area

Code
origin_pg <- st_intersection(origin_sf,pg)
origin_tm <- st_intersection(origin_sf,tm)
origin_ck <- st_intersection(origin_sf,ck)
origin_jw <- st_intersection(origin_sf,jw)

Drop redundant columns

Code
origin_pg <- origin_pg[-c(2:16)]
origin_tm <- origin_tm[-c(2:16)]
origin_ck <- origin_ck[-c(2:16)]
origin_jw <- origin_jw[-c(2:16)]

9.1.3 Preparing the lixels objects Before computing NetKDE, the SpatialLines object need to be cut into lixels with a specified minimal distance.

Code
pg_roads = st_cast(pg_roads, "LINESTRING")
tm_roads = st_cast(tm_roads, "LINESTRING")
ck_roads = st_cast(ck_roads, "LINESTRING")
jw_roads = st_cast(jw_roads, "LINESTRING")
Code
pg_lixels <- lixelize_lines(pg_roads, 
                         200, 
                         mindist = 100)
tm_lixels <- lixelize_lines(tm_roads,
                            200,
                            mindist = 100)
ck_lixels <- lixelize_lines(ck_roads,
                            200,
                            mindist = 100)
jw_lixels <- lixelize_lines(jw_roads,
                            200,
                            mindist = 100)

9.1.4 Gernerating line centre points Next, lines_center() of spNetwork will be used to generate a SpatialPointsDataFrame (i.e. samples) with line centre points

Code
pg_samples <- lines_center(pg_lixels)
tm_samples <- lines_center(tm_lixels)
ck_samples <- lines_center(ck_lixels)
jw_samples <- lines_center(jw_lixels)

9.1.5 Peforming NetKDE Puggol

Code
pg_densities <- nkde(pg_roads, 
                  events = origin_pg,
                  w = rep(1,nrow(origin_pg)),
                  samples = pg_samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5,
                  sparse = TRUE,
                  verbose = FALSE)

Tampines

Code
tm_densities <- nkde(tm_roads, 
                  events = origin_tm,
                  w = rep(1,nrow(origin_tm)),
                  samples = tm_samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5,
                  sparse = TRUE,
                  verbose = FALSE)

Choa Chu Kang

Code
ck_densities <- nkde(ck_roads, 
                  events = origin_ck,
                  w = rep(1,nrow(origin_ck)),
                  samples = ck_samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5,
                  sparse = TRUE,
                  verbose = FALSE)

Jurong West

Code
jw_densities <- nkde(jw_roads, 
                  events = origin_jw,
                  w = rep(1,nrow(origin_jw)),
                  samples = jw_samples,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5,
                  sparse = TRUE,
                  verbose = FALSE)

9.1.6 Visualising NetKDE Insert the computed density values (i.e. densities) into samples and lixels objects as density field.

Code
pg_samples$density <- pg_densities
tm_samples$density <- tm_densities
ck_samples$density <- ck_densities
jw_samples$density <- jw_densities

pg_lixels$density <- pg_densities
tm_lixels$density <- tm_densities
ck_lixels$density <- ck_densities
jw_lixels$density <- jw_densities

9.1.7 Rescalling density values

Code
pg_samples$density <- pg_samples$density*1000
tm_samples$density <- tm_samples$density*1000
ck_samples$density <- ck_samples$density*1000
jw_samples$density <- jw_samples$density*1000

pg_lixels$density <- pg_lixels$density*1000
tm_lixels$density <- tm_lixels$density*1000
ck_lixels$density <- ck_lixels$density*1000
jw_lixels$density <- jw_lixels$density*1000

Punggol

Code
tmap_mode('view')
tm_shape(pg_lixels)+
  tm_lines(col="density", lwd=5)+
  tm_shape(origin_pg)+
  tm_dots(alpha=0.5)
Code
tmap_mode('plot')

Observation : According to the Punggol map, it appears that Punggol Way has the highest density. This might be due to the limited convenience of public transportation in the area, leading residents to frequently utilize Grab services for transportation.

Tampines

Code
tmap_mode('view')
tm_shape(tm_lixels)+
  tm_lines(col="density", lwd=5)+
  tm_shape(origin_tm)+
  tm_dots(alpha=0.5)
Code
tmap_mode('plot')

Observation: According to the Tampines map, Changi Airport stands out as one of the areas with the highest density. This is logical, considering that individuals arriving in Singapore often rely on Grab services, especially since they may have luggage with them.

Choa Chu Kang

Code
tmap_mode('view')
tm_shape(ck_lixels)+
  tm_lines(col="density", lwd=5)+
  tm_shape(origin_ck)+
  tm_dots(alpha=0.5)
Code
tmap_mode('plot')

Observation : The Choa Chu Kang indicates that the central area od Choa Chu Kang, particularly in proximity to the MRT station, exhibits the highest density.

Jurong West

Code
tmap_mode('view')
tm_shape(jw_lixels)+
  tm_lines(col="density", lwd=5)+
  tm_shape(origin_jw)+
  tm_dots(alpha=0.5)
Code
tmap_mode('plot')

Observation : As per the Jurong West Map, the area with the highest density is located along Jurong West Avenue 2.